生成在C冪律分佈並用蟒

Question

測試它我知道，鑑於其產生的隨機數均勻分佈的一毫克，一種方法來獲得功率狀數據是，以下Wolfram Mathworld以下：令y是隨機可變均勻地分佈在（0,1）和x分佈爲P另一個隨機變量（x）= C * X ** N（用於（XMIN，XMAX X））。我們有

x=[ (xmax**(n+1) - xmin**(n-1))y+xmin**(n+1) ]**(1/(n+1)) 

所以我用C，其生成從1 50k的編號，以100應被分佈爲使這個節目x ^（ - 2），並打印上的文件DATA.TXT成果的頻率：

void random_powerlike(int *k, int dim, double degree, int xmin, int xmax, unsigned int *seed) 
{ 
int i; 
double aux; 
for(i=0; i<dim; i++) 
    { 
    aux=(powq(xmax, degree +1) - powq(xmin, degree +1))*((double)rand_r(seed)/RAND_MAX)+ powq(xmin, degree +1); 

    k[i]=(int) powq(aux, 1/(degree+1)); 

    } 
} 

int main() 
{ 
    unsigned int seed = 1934123471792583; 

    FILE *tmp; 
    char stringa[50]; 
    sprintf(stringa, "Data.txt"); 
    tmp=fopen(stringa, "w"); 

    int dim=50000; 
    int *k; 
    k= (int *) malloc(dim*sizeof(int)); 
    int degree=-2; 
    int freq[100]; 

    random_powerlike(k,dim, degree, 1,100,&seed); 
    fprintf(tmp, "#degree = %d x=[%d,%d]\n",degree,1,100); 
    for(int j=0; j< 100;j++) 
    { 
     freq[j]=0; 
     for(int i = 0; i< dim; ++i) 
     { 
      if(k[i]==j+1) 
      freq[j]++; 
     } 
     fprintf(tmp, "%d %d\n", j+1, freq[j]); 
    } 
    fflush(tmp); 
    fclose(tmp); 

return 0; 
} 

我決定pylab，以適應這些數字，看最好的冪律適合他們的東西作爲* X ** b，有b = -2。我在python寫了這個程序：

import numpy 
from scipy.optimize import curve_fit 
import pylab 

num, freq = pylab.loadtxt("Data.txt", unpack=True) 
freq=freq/freq[0] 

def funzione(num, a,b): 
    return a*num**(b) 

pars, covm = curve_fit(funzione, num, freq, absolute_sigma=True) 
xx=numpy.linspace(1, 99) 
pylab.plot(xx, funzione(xx, pars[0],pars[1]), color='red') 
pylab.errorbar(num, freq, linestyle='', marker='.',color='black') 
pylab.show() 
print pars 

的問題是，當我適合的數據，我得到〜-1.65的指數值。

我認爲我的地方犯了一個錯誤，但我想不出它在哪裏。

Answer 1

我認爲你必須做一個直方圖。我只是改寫你的代碼了一下，它非常適合現在

#include <math.h> 
#include <stdlib.h> 
#include <string.h> 
#include <stdio.h> 

double rndm() { 
    return (double)rand()/(double)RAND_MAX; 
} 

double power_sample(double xmin, double xmax, int degree) { 
    double pmin = pow(xmin, degree + 1); 
    double pmax = pow(xmax, degree + 1); 
    double v = pmin + (pmax - pmin)*rndm(); 
    return pow(v, 1.0/(degree + 1)); 
} 

int main() { 
    unsigned int seed = 32345U; 
    srand(seed); 

    int xmin = 1; 
    int xmax = 100; 

    double* hist = malloc((xmax-xmin + 1)*sizeof(double)); 
    memset(hist, 0, (xmax-xmin + 1)*sizeof(double)); 

    // sampling 
    int nsamples = 100000000; 
    for(int k = 0; k != nsamples; ++k) { 
     double v = power_sample(xmin, xmax, 2); 
     int idx = (int)v; 
     hist[idx] += 1.0; 
    } 

    // normalization 
    for(int k = xmin; k != xmax; ++k) { 
     hist[k] /= (double)nsamples; 
    } 

    // output 
    for(int k = xmin; k != xmax; ++k) { 
     double x = k + 0.5; 
     printf(" %e  %e\n", x, hist[k]); 
    } 

    free(hist); // cleanup 

    return 0; 
} 

及配件代碼

import numpy 
from scipy.optimize import curve_fit 
import pylab 

def funzione(x, a,b): 
    return a * numpy.power(x, b) 

num, freq = pylab.loadtxt("q.dat", unpack=True) 

pars, covm = curve_fit(funzione, num, freq, absolute_sigma=True) 
pylab.plot(num, funzione(num, pars[0], pars[1]), color='red') 
pylab.errorbar(num, freq, linestyle='', marker='.',color='black') 
pylab.show() 
print(pars) 

和它產生

[ 3.00503372e-06 1.99961571e+00] 

這是非常接近