[R適合用戶定義的分佈

我想適合我自己的分配我的數據，發現分佈相匹配的數據，最終發現在分佈中的峯的FWHM的最佳參數。從我讀過的內容來看，fitdistrplus包就是這樣做的。我知道數據在二次背景上呈現洛倫茲峯的形狀。[R適合用戶定義的分佈

圖中的數據的：使用 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)

我已經計算其限定的分佈和累積分佈方程：

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」，開始=列表（0，0.3％，-4E-04，70000，13，331））：該函數MLE未能估計參數，與在第一種情況下的錯誤代碼100'和在第二我只是得到估計值的「不適用」值列表。

我然後計算出的函數，得到位數分配值，以使用其他擬合方法（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) 
}

部分要求調用C++函數（通過使用庫RCPP）：

#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次「1」。這似乎也適用於返回明智的價值觀。分位數函數的曲線圖：

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

我選擇了「probs」值作爲隨機我不完全：

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

最後，我設法找到一個解決這個利用Python中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.

來源

2016-07-06 11:26:22 Caffrey

[R適合用戶定義的分佈

回答

相關問題