Code: The Trick


load MBHsmooths1.txt % http://www.climateaudit.info/wp-content/uploads/2009/11/mbhsmooths1.txt
% MBHsmooths1 = [ Year AnnualRecon Instrumental MB98SmoothTrick MBH98SmoothNoTrick MBH99SmoothTrick MBH99SmoothNoTrick ]

[B98,A98]=butter(10,2/50); % '50 year lowpass'

in98=MBHsmooths1(401:981,2); % Annual Recon MBH98,
%source:
% http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt ,
% col 2

Trick=MBHsmooths1(982:996,3); % pad with instrumental (1981..1995) ..
% source:
% http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt
% col 3

out98t=flipud(filter(B98,A98,flipud([filter(B98,A98,[in98;Trick;zeros(100,1)])]))); %..then smooth
% Note that filter.m initializes with zeros
out98t=out98t(1:596);
out98tr=out98t(26:576); % smoothed reconstruction with the trick

% no trick
out98nt=flipud(filter(B98,A98,flipud([filter(B98,A98,[in98;zeros(115,1)])]))); %..then smooth
out98nt=out98nt(1:596);
out98ntr=out98nt(26:576); % smoothed reconstruction without the trick

% read image, http://uc00.files.wordpress.com/2010/05/mbh985b.png
A=imread('mbh985b.png'); % from MBH98, copied for scientific purposes

XS=1.1691; % x-scale
YS=-365.7848; % y-scale
XB=-1.5749e+003; % x-bias
YB=281.8316; % y-bias
Year=(1400:1995)';

close all
image(A)

hold on
plot(Year*XS+XB,out98t*YS+YB,'r')
plot(Year*XS+XB,out98nt*YS+YB,'g')

% Note that at the begin red and green lines overlap

% MBH98: Mann, M. E., R. S. Bradley, and M. K. Hughes (1998), Global-scale
% temperature patterns and climate forcing over the past six centuries,
% Nature, 392, 779– 787.

% (C) UC 2010

2 Comments

  1. Posted May 14, 2010 at 3:14 PM | Permalink | Reply

    You might just have to open that up to people. Just the source files and and the results with anything notable to see would be great.

    If there were a single line at the start, stating “this is the code that fooled the world” might make a difference, but that’s something that one won’t be looking for. The problem, as it has always been, is usually a touch more subtle.

    -Sale

  2. Steve McIntyre
    Posted Mar 16, 2011 at 1:32 PM | Permalink | Reply

    Here is implementation of UC’s code in R:

    ##LOAD MBH data
    #load MBHsmooths1.txt # http://www.climateaudit.info/wp-content/uploads/2009/11/mbhsmooths1.txt
    ## MBHsmooths1 = [ Year AnnualRecon Instrumental MB98SmoothTrick MBH98SmoothNoTrick MBH99SmoothTrick MBH99SmoothNoTrick ]

    MBH=MBHsmooths1=read.table(“http://www.climateaudit.info/wp-content/uploads/2009/11/mbhsmooths1.txt”)
    names(MBH)=names(MBHsmooths1)=c(“Year”, “AnnualRecon”, “Instrumental” ,”MBH98SmoothTrick”, “MBH98SmoothNoTrick”,
    “MBH99SmoothTrick”, “MBH99SmoothNoTrick”)
    dim(MBH) #999 7
    MBH[1,]

    #[B98,A98]=butter(10,2/50); # ’50 year lowpass’

    library(signal)
    bf=butter(10,2/50); # ’50 year lowpass’
    B98=bf$b;A98=bf$a
    B98 # [1] 6.500312e-13 6.500312e-12 2.925140e-11

    #in98=MBHsmooths1(401:981,2); # Annual Recon MBH98,
    ##source:
    ## http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt ,
    ## col 2

    in98= MBHsmooths1[401:981,"AnnualRecon"];
    length(in98) #581

    #Trick=MBHsmooths1(982:996,3) #; # pad with instrumental (1981..1995) ..
    ## source:
    ## http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt
    ## col 3

    Trick=MBHsmooths1[982:996,"Instrumental"]#; # pad with instrumental (1981..1995) ..
    length(Trick) #15

    ##function flipud
    #tfilt=conv(E,flipud(E))/M;
    #Oh, I forgot to explain the flipud-thing:
    #Suppose E=[1;2;3;4;5] (a column vector)
    # Then flipud simply reverses the values (each column), i.e. flipud(E)=[5;4;3;2;1] and

    flipud=function(x) rev(x)

    # conv(E,flipud(E)) is
    # 5 14 26 40 55 40 26 14 5

    E=1:5
    conv(E,flipud(E))

    #They should be convolved
    #tfilt=conv(E,flipud(E))/M;
    #Reading the point 145 in the document I referred that would be in R something like tfilt=convolve(E,E,type=’open’)/M

    convolve(E,E,type=”open”)
    # 5 14 26 40 55 40 26 14 5

    #out98t=flipud(filter(B98,A98,flipud( [filter(B98,A98, [in98;Trick;zeros(100,1)] )] ))); #..then smooth
    #% Note that filter.m initializes with zeros
    #out98t=out98t(1:596);
    #out98tr=out98t(26:576); % smoothed reconstruction with the trick
    #flipud(filter(B98,A98,flipud( [filter(B98,A98, c(in98,Trick, rep(0,100) ) )] )));
    x= c( in98,Trick )
    #flipud(filter(B98,A98,flipud( [ filter(B98,A98, c(in98,Trick, rep(0,100) ) ) ] )));
    # y= filter(B98,A98, x )
    out98t= y= rev (filter(B98,A98,rev ( filter(B98,A98,c(rep(0,100),x,rep(0,100) ) ) ) ))

    out98t=flipud(filter(B98,A98, flipud( c(filter(B98,A98, c(in98,Trick,rep(0,100)) )) ))); #..then smooth
    length(out98t)
    out98t=ts(out98t[1:596],start=1400)
    out98tr=ts( out98t[26:576],start=1425) # % smoothed reconstruction with the trick

    plot(1400:1980,in98,type=”l”,col=”grey50″,ylim=c(-.8,.5))
    lines(out98tr,lty=2,lwd=3)

    cbind(out98tr ,MBH[(1425:1975)-999,"MBH98SmoothTrick"]) #1425 to 1975
    range( out98tr -MBH[(1425:1975)-999,"MBH98SmoothTrick"],na.rm=TRUE)
    # -0.000729800 0.001258341

    # no trick
    #out98nt=flipud(filter(B98,A98,flipud([filter(B98,A98,[in98;zeros(115,1)])]))); %..then smooth
    #out98nt=out98nt(1:596);
    #out98ntr=out98nt(26:576); % smoothed reconstruction without the trick

    out98nt=flipud(filter(B98,A98, flipud( c(filter(B98,A98, c(in98,rep(0,115)) )) ))); #..then smooth
    out98nt=ts(out98nt[1:596],start=1400)
    out98ntr=ts( out98nt[26:576],start=1425) # % smoothed reconstruction with the trick

    range( out98ntr -MBH[(1425:1975)-999,"MBH98SmoothNoTrick"],na.rm=TRUE)
    # -5.796738e-05 1.489949e-04

    #% read image, http://uc00.files.wordpress.com/2010/05/mbh985b.png
    #A=imread(‘mbh985b.png’); % from MBH98, copied for scientific purposes
    #XS=1.1691; % x-scale
    #YS=-365.7848; % y-scale
    #XB=-1.5749e+003; % x-bias
    #YB=281.8316; % y-bias
    #Year=(1400:1995)';
    #close all
    #image(A)

    library(png)
    download.file(“http://uc00.files.wordpress.com/2010/05/mbh985b.png”,”temp.png”,mode=”wb”)
    img=readPNG(“temp.png”)
    dim(img)
    #610 781 3

    XS=1.1691; # % x-scale
    YS=-365.7848; #% y-scale
    XB=-1.5749e+003;# % x-bias
    YB=281.8316;# % y-bias
    Year=(1400:1995);

    library(grDevices)
    par(mar=c(0,0,0,0))
    par(bg=”lightblue”)
    plot(MBH[,1],xlim=c(1300,2025),ylim=c(-.8,.5),type=”n”,xaxs=”i”,axes=FALSE)
    rasterImage(img , 1320, -.65, 2025, .5)
    lines(1275+(Year-1300)*1.05,out98nt*.7-.03,col=3,lwd=1)
    lines(1275+(Year-1300)*1.05,out98t*.7-.03,col=2,lwd=1)

    #temp=Year>=1980
    #lines(1275+(Year[temp]-1300)*1.05,out98t[temp]*.7-.03,col=2,lwd=2)
    temp=Year>=1940
    lines(1275+(Year[temp]-1300)*1.05,out98nt[temp]*.7-.03,col=”yellow”,lwd=2)

3 Trackbacks

  1. By The Trick Timeline « Climate Audit on Jul 1, 2010 at 11:21 AM

    [...] Fools, here’s the turn-key(*) [...]

  2. [...] diagnosed by CA reader UC here and expounded in greater length (with Matlab code here and here and here ). It consists of the following elements: 1. A digital splice of proxy data up to 1980 [...]

  3. [...] diagnosed by CA reader UC here and expounded in greater length (with Matlab code here and here and here ). It consists of the following elements: 1. A digital splice of proxy data up to 1980 [...]

Post a Comment

Required fields are marked *

*
*

Follow

Get every new post delivered to your Inbox.

Join 3,136 other followers

%d bloggers like this: