如何按自定义的坐标点生成随机抽样

目录

1.问题描述

2.解决方法

3.代码实现

1.问题描述

更具自定义的数据点，生成能量的随机抽样；

需要的数据分布点；

2.解决方法

先积分然后归一化得到概率函数，采用二分法找到对应的区间的X；

生成的随机抽样分布如上；

3. 代码实现

#include "iostream"
using namespace std;

void GenerateTest()
{
  vector<double> bins,vals,nvals; 
  //读取数据
  double x=0.,y=0.;
  ifstream input("Alt100.txt");

  bins.clear();vals.clear();
  bins.push_back(0);
  vals.push_back(0);
  for(int i=1;i<191;i++)
    {
      input>>x>>y;
      bins.push_back(x);
      vals.push_back(y);
    }

  //求和和归一化
  double sum=0;
  vals.clear();
  for(int ii=1;ii<191;ii++)
    {
      // value=vals[ii]*(bins[ii]-bins[ii-1]);
      // sum = sum+value;
      vals[ii]=vals[ii];//*(bins[ii]-bins[ii-1]);
      sum = sum+vals[ii];
      // cout<<"  Sum"<<sum<<endl;
      // cout<<""<<vals[ii]<<endl;
    }

  // for(int ii=1;ii<190;ii++)
  //   {
  //     vals[ii]=vals[ii]/(bins[ii]-bins[ii-1]);
  //     // sum = sum+vals[ii];
  //     // cout<<"  Sum"<<sum<<endl;
  //     // cout<<""<<vals[ii]<<endl;
  //   }

  nvals.clear();
  // 
  double norm=0.;
  for(int iii=1;iii<191;iii++)
    {
      norm = vals[iii]/sum+norm;
      cout<<"  norm "<<norm<<endl;
      // cout<<" vals"<<vals[iii]<<endl;
      nvals.push_back(norm);
    }

  // for(int j=0;j<189;j++) cout<<nvals[j]<<endl;//test nvals
  //
  // double Emin=1.0e-8;
  // double Emax=1000;

  TRandom3 rndm;
  //--------对数坐标变化--------
  double low=1.e-8;
  double up=1e4;
  double nbin=192;
  Double_t xAxis[193];

  for ( int i=0; i <= nbin; i++) {
    double val_bin=low * std::pow(10., i * std::log10(up/low)/nbin);
    double exp_10=4.-int(std::log10(val_bin));
    double factor =std::pow(10., exp_10);
    val_bin=int(factor*val_bin)/factor;
    xAxis[i] = val_bin;
  }


  TH1D* h1 = new TH1D("h1","h1",192,xAxis);
  // TH1D* h1 = new TH1D("h1","h1",400,0.0,1.0);
  //--------------------------
  for(int j=0;j<1000000;j++){
    int numberOfBin=nvals.size();
    int nmin=0;
    int nmid=numberOfBin/2;
    int nmax=numberOfBin-1;
    double a=rndm.Uniform();
    while(nmin!=nmax-1)
      {
	if(a>nvals[nmid])
	  nmin=nmid;
	else
	  nmax=nmid;
	nmid=nmin+(nmax-nmin+1)/2;
      }

    double weight=(bins[nmid]-bins[nmid-1])/(nvals[nmid]-nvals[nmid-1]);
    double rand=bins[nmid-1]+weight*(a-nvals[nmid-1]);
    h1->Fill(rand);

  }

  TCanvas *c1 = new TCanvas();
  c1->cd();
  TGraph *gr=new TGraph("./Alt100.txt","%lg%lg");
  gr->Draw();
  TCanvas *c2 = new TCanvas();
  c2->cd();
  //
  h1->Draw();
}