import pandas as pdimport numpy as np= np.random.default_rng(seed= 5650 )import matplotlib.pyplot as pltimport seaborn as sns= ['#e69f00' ,'#56b4e9' ,'#009e73' ]import patchworklib as pwfrom scipy.special import expit# The original R code: # sim_happiness <- function( seed=1977 , N_years=1000 , max_age=65 , N_births=20 , aom=18 ) { # set.seed(seed) # H <- M <- A <- c() # for ( t in 1:N_years ) { # A <- A + 1 # age existing individuals # A <- c( A , rep(1,N_births) ) # newborns # H <- c( H , seq(from=-2,to=2,length.out=N_births) ) # sim happiness trait - never changes # M <- c( M , rep(0,N_births) ) # not yet married # # for each person over 17, chance get married # for ( i in 1:length(A) ) { # if ( A[i] >= aom & M[i]==0 ) { # M[i] <- rbern(1,inv_logit(H[i]-4)) # } # } # # mortality # deaths <- which( A > max_age ) # if ( length(deaths)>0 ) { # A <- A[ -deaths ] # H <- H[ -deaths ] # M <- M[ -deaths ] # } # } # d <- data.frame(age=A,married=M,happiness=H) # return(d) # DGP: happiness -> marriage <- age = 70 = 41 = ['age' ,'a' ,'h' ,'m' ]= []= np.zeros(shape= (num_births,1 ))= np.linspace(- 2 , 2 , num= num_births)= np.zeros(shape= (num_births,1 ))def update_m(row):if row['m' ] == 0 :return int (rng.binomial(= 1 ,= expit(row['h' ] - 3.875 ),= 1 ,return 1 def sim_cohort_to(max_age):= pd.DataFrame({'age' : [1 for _ in range (num_births)],'h' : np.linspace(- 2 , 2 , num= num_births),'m' : [0 for _ in range (num_births)],for t in range (2 , max_age + 1 ):'age' ] = sim_df['age' ] + 1 if t >= 18 :'m' ] = sim_df.apply (update_m, axis= 1 )return sim_df= []for cur_max_age in range (1 , 71 ):= sim_cohort_to(cur_max_age)= pd.concat(all_sim_dfs)# full_sim_df.head() = ['#c6c6c666' ] + cb_palette'm_label' ] = full_sim_df['m' ].apply (lambda x: "Unmarried" if x == 0 else "Married" )= full_sim_df.rename(columns= {'age' : 'Age' , 'h' : 'Happiness' })= pw.Brick(figsize= (5.25 ,2.75 )); = 'Age' , y= 'Happiness' , hue= 'm_label' ,= full_sim_df,= ax,= cbg_palette,= 22 ,= True ,; "upper center" , bbox_to_anchor= (0.5 , 1.15 ), ncol= 2 ); "" ); = 17.5 , color= 'black' , ls= 'dashed' , lw= 1 ); = full_sim_df.groupby('Age' )['Happiness' ].mean().reset_index()'m_label' ] = "All" 'Happiness_mean' ] = 0 = full_sim_df[full_sim_df['Age' ] >= 18 ].groupby(['Age' ,'m_label' ])['Happiness' ].mean().reset_index()= group_hap_df[group_hap_df['m_label' ] == "Married" ].copy()= group_hap_df[group_hap_df['m_label' ] == "Unmarried" ].copy()# Moving averages = 3 'Happiness_mean' ] = married_df['Happiness' ].rolling(window= win_size).mean()'Happiness_mean' ] = unmarried_df['Happiness' ].rolling(window= win_size).mean()# Merge back together = pd.concat([mean_hap_df, married_df, unmarried_df], ignore_index= True )# display(combined_df.tail()) = pw.Brick(figsize= (3.75 ,2.4 )); 0 , 70 ); # ax.set_ylim(-1.2, 1.2); = ['black' , cb_palette[0 ], '#b2b2b2ff' ]# And plot = 'Age' , y= 'Happiness_mean' ,= 'm_label' ,# marker=".", # s=20, = combined_df,= cbg_palette,# color=cbg_palette[0], # order=3, = ax; "Mean Happiness by Status" ); "" ); "Mean Happiness" ); 
