生成在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] 

这是非常接近