How can I manually calculate bpy.types.PoseBone.matrix using Blender's Python API?
armature.data.bones["Bone"].matrix_local


Is the bone's world un-posed matrix (the pose_bone.matrix in rest state).

 armature.pose.bones["Bone"].matrix_basis


Is the bone's local pose matrix.

From there we can calculate any bone's world matrix by recursively traversing it's parents. Note: This will not include any constraints or IK (but it will include drivers or actions):

def matrix_world(armature_ob, bone_name):
    local = armature_ob.data.bones[bone_name].matrix_local
    basis = armature_ob.pose.bones[bone_name].matrix_basis

    parent = armature_ob.pose.bones[bone_name].parent
    if parent == None:
        return  local * basis
    else:
        parent_local = armature_ob.data.bones[parent.name].matrix_local
        return matrix_world(armature_ob, parent.name) * (parent_local.inverted() * local) * basis


You can test it on an empty and some armature with some active pose_bone. Also the armature_object's transformation is accounted here:

import bpy

empty_ob = bpy.data.objects["Empty"]
armature_ob = bpy.data.objects["Armature"]

bone_name = bpy.context.active_pose_bone.name
empty_ob.matrix_world = armature_ob.matrix_world * matrix_world(armature_ob, bone_name )