Generating a power-law distribution in C and testing it with python -


i know that, given rng generates random numbers uniformly distributed, way obtain power-like data is, following wolfram mathworld following: let y random variable uniformly distributed in (0,1) , x random variable distributed p(x) = c*x**n (for x in (xmin,xmax)). have that

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

so made program in c generates 50k numbers 1 100 should distributed x^(-2) , prints frequency of outcomes on file 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 = 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; }

i decided fit these numbers pylab, see if best power-law fit them a*x**b, b = -2. wrote program in python:

import numpy 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

the problem when fit data, obtain exponent value of ~-1.65.

i think made mistake somewhere, can't figure out where.

i think have make histogram. rewrote code bit , fits now

#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; }

and fitting code

import numpy scipy.optimize import curve_fit import pylab  def funzione(x, a,b):     return * 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)

and produced 

[  3.00503372e-06   1.99961571e+00]

which pretty close


Comments

Post a Comment

Popular posts from this blog

android - InAppBilling registering BroadcastReceiver in AndroidManifest -

i'm implementing inapp billing in android app. works well, however, i'm trying decouple broadcast receiver activity manifest. suggestion in android's trivialdrive sample: // important: dynamically register broadcast messages updated purchases. // register receiver here instead of <receiver> in manifest // because call getpurchases() @ startup, therefore can ignore // broadcasts sent while app isn't running. // note: registering listener in activity bad idea, done here // because sample. regardless, receiver must registered after // iabhelper setup, before first call getpurchases(). currently there's class extends broadcastreceiver : public class iabbroadcastreceiver extends broadcastreceiver { /** * intent action receiver should filter for. */ public static final string action = "com.android.vending.billing.purchases_updated"; private final iabbroadcastlistener mlistener; public iabbroadcastreceiver(iabbroadcastlistener listener) { m...
Read more

python Tkinter Capturing keyboard events save as one single string -

is possible store keyboard events 1 string? code below stores 1 char , prints it. card reader or bar code reader, contains collection of character/ string not 1 character @ time. goal save char pressed text variable. from tkinter import * root = tk() def key(event): text= event.char text+= event.char print ("pressed", text) def callback(event): frame.focus_set() print ("clicked at", event.x, event.y) frame = frame(root, width=100, height=100) frame.bind("<key>", key) frame.bind("<button-1>", callback) frame.pack() root.mainloop() currently, creating text variable , everytime key() function gets called, , text stores last character typed. you can define text module level variable , , use module level text inside key function - text = '' def key(event): global text text+= event.char print("pressed", text)
Read more

sql server - Why does Linq-to-SQL add unnecessary COUNT()? -

here tables: create table customer ( cid int primary key, name varchar(32) ) create table ord ( oid int primary key, cid int foreign key references customer, address varchar(20) ) here's linq-to-sql statement: var aaa = c in db.customer select new { c, o = c.ord.tolist() }; and here's inexplicable query generated linq-to-sql: select [t0].[cid] [cid], [t0].[name] [name], [t1].[oid] [oid], [t1].[cid] [cid2], [t1].[address] [address], (select count(*) [dbo].[ord] [t2] [t2].[cid] = [t0].[cid]) [value] [dbo].[customer] [t0] left outer join [dbo].[ord] [t1] on [t1].[cid] = [t0].[cid] order [t0].[cid], [t1].[oid] i want understand how rid of count(*) part. it's entirely uncalled for! linq needs count(*) in order determine number of entries returned o = c.ord.tolist() .
Read more