Calculating final Bone/PoseBone matrices
A straight-forward approach using convert_space() and a matrix to convert bone space:

import bpy

arma_ob = bpy.context.object
mesh_ob = bpy.data.objects['Suzanne']

pose_bone = bpy.context.active_pose_bone

bone_space = Matrix(((1,0,0,0),(0,0,1,0),(0,-1,0,0),(0,0,0,1)))
mat = arma_ob.convert_space(pose_bone, pose_bone.matrix, 'POSE', 'WORLD') * bone_space

mesh_ob.matrix_world = mat


Add a monkey and select a pose bone, then run the above script. Then Suzanne should look in pose bone direction. Note that bone direction is +Y, so it's (0,1,0) for a bone pointing upwards (Z).