Rフィットユーザー定義分布

Question

データに自分自身の分布を適合させようとしています。データと一致する分布の最適なパラメーターを見つけ、最終的には分布のピークの半値全幅を求めます。私が読んだことから、fitdistrplusというパッケージがこれを行う方法です。私は、データが2次バックグラウンド上のローレンツピークの形をとっていることを知っています。データの

プロット：使用 plot of raw data

生データ：

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は、分布および累積分布を定義する式を計算した：

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)) 

私はこれらを信じ正しいこと。私が理解から、分布フィットfitdist（またはmledist）を使用して、ちょうどこれらの定義を使用可能にする必要があります方法：

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)) 

これは、文「関数を返す初期パラメータで評価することはできません> fitdistでエラーが発生しました（データ、 "FF"、start = list（0,0.3、-4e-04,70000,13,331））：関数mleはパラメータの推定に失敗し、エラーコード100'が第1の場合と第2の場合私はちょうど推定のための「NA」の値のリストを得ます。

Iは、他のフィッティング方法（qmefit）を使用する分位数分布値を与える関数を計算：このコードの

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) 
} 

一部（ライブラリRcppを使用して）、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の場合、返されたベクトルなどに4が加算されます。これはまた、合理的な値が返されるときに機能するようです。クォン関数のプロット：私は完全にはないそうであるようにランダムに「ちゃったごめんなさい」の値を選択しました

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>

私はここで基本的な間違いを犯しているかどうかわかりません。感謝。

Answer 1

最後に、私はrPythonパッケージとlmfitを使ってこれを回避する方法を見つけました。それは私の問題を解決し、同じ問題を抱える他の人にも役立つかもしれません。次のようにR-コードはでした：

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") 

と呼ばれるPythonのコードはでした：

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 

私はその誤り、データのためのより良いフィット感を得るために、多項式の次数を増加させ、FWHMを返さフィットの値。これを達成するにはもっと良い方法がありますが、最終的には必要なものがあります。

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