R 적합 사용자 정의 분포

데이터에 내 자신의 분포를 맞추려고합니다. 데이터와 일치하는 최적의 분포 매개 변수를 찾고 분포에서 최고점의 반 누름을 찾습니다. 필자가 읽은 바로는 fitdistrplus 패키지가이를 수행하는 방법입니다. 나는 데이터가 2 차 배경에 lorentzian peak의 형태를 취하고 있음을 안다. 데이터R 적합 사용자 정의 분포

플롯 : plot of raw data

원시 데이터를 사용 : 이러한 믿을

dFF <- function(x,a,b,c,A,gamma,pos) a + b*x + (c*x^2) + ((A/pi)*(gamma/(((x-pos)^2) + (gamma^2)))) 

pFF <- function(x,a,b,c,A,gamma,pos) a*x + (b/2)*(x^2) + (c/3)*(x^3) + A/2 + (A/pi)*(atan((x - pos)/gamma))

data = c(0,2,5,4,5,4,3,3,2,2,0,4,4,2,5,5,3,3,4,4,4,3,3,5,5,6,6,8,4,0,6,5,7,5,6,3,2,1,7,0,7,9,5,7,5,3,5,5,4,1,4,8,10,2,5,8,7,14,7,5,8,4,2,2,6,5,4,6,5,7,5,4,8,5,4,8,11,9,4,8,11,7,8,6,9,5,8,9,10,8,4,5,8,10,9,12,10,10,5,5,9,9,11,19,17,9,17,10,17,18,11,14,15,12,11,14,12,10,10,8,7,13,14,17,18,16,13,16,14,17,20,15,12,15,16,18,24,23,20,17,21,20,20,23,20,15,20,28,27,26,20,17,19,27,21,28,32,29,20,19,24,19,19,22,27,28,23,37,41,42,34,37,29,28,28,27,38,32,37,33,23,29,55,51,41,50,44,46,53,63,49,50,47,54,54,43,45,58,54,55,67,52,57,67,69,62,62,65,56,72,75,88,87,77,70,71,84,85,81,84,75,78,80,82,107,102,98,82,93,98,90,94,118,107,113,103,99,103,96,108,114,136,126,126,124,130,126,113,120,107,107,106,107,136,143,135,151,132,117,118,108,120,145,140,122,135,153,157,133,130,128,109,106,122,133,132,150,156,158,150,137,147,150,146,144,144,149,171,185,200,194,204,211,229,225,235,228,246,249,238,214,228,250,275,311,323,327,341,368,381,395,449,474,505,529,585,638,720,794,896,919,1008,1053,1156,1134,1174,1191,1202,1178,1236,1200,1130,1094,1081,1009,949,890,810,760,690,631,592,561,515,501,489,467,439,388,377,348,345,310,298,279,253,257,259,247,237,223,227,217,210,213,197,197,192,195,198,201,202,211,193,203,198,202,174,164,162,173,170,184,170,168,175,170,170,168,162,149,139,145,151,144,152,155,170,156,149,147,158,171,163,146,151,150,147,137,123,127,136,149,147,124,137,133,129,130,128,139,137,147,141,123,112,136,147,126,117,116,100,110,120,105,91,100,100,105,92,88,78,95,75,75,82,82,80,83,83,66,73,80,76,69,81,93,79,71,80,90,72,72,63,57,53,62,65,49,51,57,73,54,56,78,65,52,58,49,47,56,46,43,50,43,40,39,36,45,28,35,36,43,48,37,36,35,39,31,24,29,37,26,22,36,33,24,31,31,20,30,28,23,21,27,26,29,21,20,22,18,19,19,20,21,20,25,18,12,18,20,20,13,14,21,20,16,18,12,17,20,24,21,20,18,11,17,12,5,11,13,16,13,13,12,12,9,15,13,15,11,12,11,8,13,16,16,16,14,8,8,10,11,11,17,15,15,9,9,13,12,3,11,14,11,14,13,8,7,7,15,12,8,12,14,9,5,2,10,8)

는 I의 분포와 누적 분포를 정의하는 방정식을 산출했다 올바른지. 내가 이해에서 배포 적합이 fitdist (또는 mledist) 메소드를 사용하여 바로 이러한 정의 사용 가능해야한다 : 이것은 fitdist에서 '기능 초기 매개 변수에서 평가 될 수없는 문> 오류를 반환

fitdist(data,'FF', start = list(0,0.3,-0.0004,70000,13,331)) mledist(data,'FF', start = list(0,0.3,-0.0004,70000,13,331))

을 (데이터 , "FF", 시작 = 목록 (0, 0.3, -4e-04, 70000, 13, 331)) : 함수 mle은 매개 변수를 추정하지 못했습니다. 첫 번째 경우와 두 번째 경우에 ' 나는 견적을위한 'NA'값의 목록을 얻는다.

그때 다른 결합 방법 (qmefit) 사용하는 분위수 분포 값을 부여하는 함수를 계산

(라이브러리 Rcpp 사용)이 코드의

qFF <- function(p,a,b,c,A,gamma,pos) 
{ 
    qList = c() 
    axis = seq(1,600,1) 
    aF = dFF(axis,a,b,c,A,gamma,pos) 
    arr = histogramCpp(aF) # change data to a histogram format 
    for(element in 1:length(p)){ 
     q = quantile(arr,p[element], names=FALSE) 
     qList = c(qList,q) 
    } 
    return(qList) 
}

부분은 C++ 함수를 호출 할 필요가 :

#include <Rcpp.h> 
#include <vector> 
#include <math.h> 
using namespace Rcpp; 

// [[Rcpp::export]] 
std::vector<int> histogramCpp(NumericVector x) { 
    std::vector<int> arr; 
    double number, fractpart, intpart; 
    for(int i = 0; i <= 600; i++){ 
     number = (x[i]); 
     fractpart = modf(number , &intpart); 
     if(fractpart < 0.5){ 
      number = (int) intpart; 
     } 
     if(fractpart >= 0.5){ 
      number = (int) (intpart+1); 
     } 
     for(int j = 1; j <= number; j++){ 
      arr.push_back(i); 
     } 
    } 
    return arr; 
}

이 C++ 메소드는 데이터를 히스토그램 형식으로 변환합니다. 데이터를 설명하는 벡터의 첫 번째 요소가 4이면 '1'이 반환 된 벡터 등에 4 번 추가됩니다. 이것은 또한 합리적인 값이 반환되기 때문에 효과가있는 것으로 보입니다. 분위수 기능의 줄거리 : 나는 완벽하지 않기 때문에 무작위로 'probs'값을 선택

fitdist(data,'FF', start = list(0,0.3,-0.0004,70000,13,331), method = 'qme', probs = c(0,0.3,0.4,0.5,0.7,0.9))

Plot of quantiles returned for probabilities from 0 to 1 in steps of 0.001

'qmefit'방법은 다음 fitdist 기능을 통해 시도 할 수 있습니다 그들의 의미를 이해하십시오. 이것은 곧장 R 세션이 충돌하거나 간단한 말더듬이 'NA'값의리스트를 견적과 회신으로 반환하면 <std::bad_alloc : std::bad_alloc>

내가 여기에 기본적인 실수를하고 있는지, 어떤 도움이나 권고 사항이 고맙습니다.

출처

2016-06-29 Caffrey

결국 나는 rPython 패키지와 lmfit에서 파이썬을 사용하여이 문제를 해결할 수있었습니다. 그것은 내 문제를 해결하고 같은 문제가있는 다른 사람들에게 유용 할 수 있습니다.

library(rPython) 
python.load("pyFit.py") 
python.assign("row",pos) 
python.assign("vals",vals) 
python.exec("FWHM,ERROR,FIT = fitDist(row,vals)") 
FWHM = python.get("FWHM") 
ERROR = python.get("ERROR") 
cFIT = python.get("FIT")

와 호출 된 파이썬 코드가 있었다 : 나는 그 오류 데이터에 대한 더 적합을 얻기 다항식의 순서를 증가하고 FWHM을 반환

from lmfit import Model, minimize, Parameters, fit_report 
from sklearn import mixture 
import numpy as np 
import matplotlib.pyplot as plt 
import math 

def cauchyDist(x,a,b,c,d,e,f,g,A,gamma,pos): 
    return a + b*x + c*pow(x,2) + d*pow(x,3) + e*pow(x,4) + f*pow(x,5) + g*pow(x,6) + (A/np.pi)*(gamma/((pow((x-pos),2)) + (pow(gamma,2)))) 

def fitDist(row, vals): 
    gmod = Model(cauchyDist) 

    x = np.arange(0,600) 
    result = gmod.fit(vals, x=x, a = 0, b = 0.3, c = -0.0004, d = 0, e = 0, f= 0, g = 0, A = 70000, gamma = 13, pos = row) 

    newFile = open('fitData.txt', 'w') 
    newFile.write(result.fit_report()) 
    newFile.close() 

    with open('fitData.txt', 'r') as inF: 
     for line in inF: 
      if 'gamma:' in line: 
       j = line.split() 
    inF.close() 

    FWHM = float(j[1]) 
    error = float(j[3]) 
    fit = result.best_fit 
    fit = fit.tolist() 
    return FWHM, error, fit

다음과 같이 R-코드이었다 그리고 적합성에 대한 값. 이것을 달성하는 훨씬 더 좋은 방법이있을 수 있지만 최종 적합은 내가 필요로하는 것입니다.

Final fit. Red data points are raw data, the black line is the fitted distribution.

출처

2016-07-06 11:26:22 Caffrey

답변

관련 문제