Ignore:
Timestamp:
Oct 7, 2013 12:30:10 PM (8 years ago)
Author:
vondreele
Message:

simulated CW powder data now complete with Poisson errors

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r1078 r1087  
    1818import numpy.linalg as nl
    1919import scipy.optimize as so
     20import scipy.stats as st
    2021import GSASIIpath
    2122GSASIIpath.SetVersionNumber("$Revision$")
     
    18761877            Limits = calcControls[hfx+'Limits']
    18771878            x,y,w,yc,yb,yd = Histogram['Data']
    1878             W = wtFactor*w
    18791879            yc *= 0.0                           #zero full calcd profiles
    18801880            yb *= 0.0
     
    18821882            xB = np.searchsorted(x,Limits[0])
    18831883            xF = np.searchsorted(x,Limits[1])
     1884            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
     1885                varylist,Histogram,Phases,calcControls,pawleyLookup)
     1886            yc[xB:xF] += yb[xB:xF]
     1887            if not np.any(y):                   #fill dummy data
     1888                rv = st.poisson(yc[xB:xF])
     1889                y[xB:xF] = rv.rvs()
     1890                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
     1891            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
     1892            W = wtFactor*w
     1893            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
    18841894            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
    18851895            Nobs += Histogram['Residuals']['Nobs']
    18861896            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
    18871897            SumwYo += Histogram['Residuals']['sumwYo']
    1888             yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
    1889                 varylist,Histogram,Phases,calcControls,pawleyLookup)
    1890             yc[xB:xF] += yb[xB:xF]
    1891             yd[xB:xF] = y[xB:xF]-yc[xB:xF]
    1892             wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
    18931898            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
    18941899            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
Note: See TracChangeset for help on using the changeset viewer.